/*
 * Reaction.cpp
 *
 *  Created on: 23 May 2011
 *      Author: allan
 */

#include "Reaction.h"

// C++ includes
#include <cmath>
#include <iomanip>

// Boost includes
#include <boost/foreach.hpp>

// GeoReact includes
#include "NumericalUtilities.h"

ReactionEquation::ReactionEquation(const SpeciesMap& speciesMap) :
speciesMap(speciesMap)
{}

void ReactionEquation::Set(const string& equation)
{
	// TODO
}

void ReactionEquation::Set(const string& reactant, const vector<string>& products, const vector<double>& stoichiometries)
{
	// Set the reactant species of the reaction
	this->reactant = {reactant, speciesMap[reactant], 1.0};
	
	// Set the product species of the reaction
	for(unsigned i = 0; i < products.size(); ++i)
	{
		// Create a new product species
		ReactionSpecies product = {products[i], speciesMap[products[i]], stoichiometries[i]}; 

		// Insert the new product into (products)
		this->products.push_back(product);
	}
}

const ReactionSpecies& ReactionEquation::GetReactant() const
{
	return reactant;
}

const vector<ReactionSpecies>& ReactionEquation::GetProducts() const
{
	return products;
}

Reaction::Reaction()
{}

Reaction::Reaction(const ReactionEquation& equation) :
reactant(equation.GetReactant()), products(equation.GetProducts())
{
	SetKinectics(ZeroKinetics());
}
	
Reaction::Reaction(const ReactionEquation& equation, const EquilibriumConstant& K) :
reactant(equation.GetReactant()), products(equation.GetProducts()), K(K)
{
	SetKinectics(ZeroKinetics());
}

void Reaction::SetEquilibriumConstant(const EquilibriumConstant& K)
{
	this->K = K;
}

const ReactionSpecies& Reaction::GetReactant() const
{
	return reactant;
}

const vector<ReactionSpecies>& Reaction::GetProducts() const
{
	return products;
}

const EquilibriumConstant& Reaction::GetEquilibriumConstant() const
{
	return K;
}

const double Reaction::ReactionQuotient(const VectorXd& a) const
{
	// The reaction quotient
	double Qr = 1.0/a[reactant.index];
	
	for(unsigned i = 0; i < products.size(); ++i)
	{
		const double ai = a[products[i].index];
		const double vi = products[i].stoichiometry;
		
		Qr *= std::pow(ai, vi);
	}
	
	return Qr;
}

const double Reaction::MassRate(double T, double P, const VectorXd& n, const VectorXd& a) const
{
	return kinetics->MassRate(T, P, n, a);
}

ostream& operator<<(ostream& out, const Reaction& reaction)
{
	out << setw(10) << reaction.GetReactant().name << " (" << reaction.GetReactant().index << ")" << " = ";
	
	BOOST_FOREACH(auto product, reaction.GetProducts())
		out << setprecision(2) << setw(4) << product.stoichiometry << setw(10) << product.name << " (" << product.index << ")" << ", ";
	
	return out;
}

